/*******************************************************************************
* pymesh. mesh routines for nemesis finite element code.                       *
* Copyright (C) 2009-2011 F.E.Karaoulanis [http://www.nemesis-project.org]     *
*                                                                              *
* This program is free software; you can redistribute it and/or modify         *
* it under the terms of the GNU General Public License version 3, as           *
* published by the Free Software Foundation.                                   *
*                                                                              *
* This program is distributed in the hope that it will be useful,              *
* but WITHOUT ANY WARRANTY; without even the implied warranty of               *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                *
* GNU General Public License for more details.                                 *
*                                                                              *
* You should have received a copy of the GNU General Public License            *
* along with this program.  If not, see < http://www.gnu.org/licenses/>.       *
*******************************************************************************/

// *****************************************************************************
// $LastChangedDate$
// $LastChangedRevision$
// $LastChangedBy$
// $HeadURL$
// Author(s): F.E. Karaoulanis (fkar@nemesis-project.org)
// *****************************************************************************

#ifdef SINGLE
#define REAL float
#else /* not SINGLE */
#define REAL double
#endif /* not SINGLE */

#include <Python.h>
#include <structmember.h>

#define REAL double
#define VOID int
#include "triangle/triangle.h"
#define DEBUG 1
#define DEBUG_PRINT(fmt, ...) \
        do { if (DEBUG) fprintf(stderr, "%s:%d:%s(): " fmt, __FILE__, \
                                __LINE__, __func__, __VA_ARGS__); } while (0)

/*******************************************************************************
 * Data structure definition and functions
 ******************************************************************************/
/**
 * Data structure.
 */
typedef struct {
  PyObject_HEAD
  /* Type-specific fields go here. */
  struct triangulateio* io;
  int pointlistsize;
  int pointattributelistsize;
  int pointmarkerlistsize;
  int edgelistsize;
  int edgemarkerlistsize;
  int holelistsize;
  int neighborlistsize;
  int normlistsize;
  int regionlistsize;
  int segmentmarkerlistsize;
  int segmentlistsize;
  int trianglelistsize;
  int triangleattributelistsize;
  int trianglearealistsize;
 } MeshIO;

/**
 * Deallocate All memory used. Initialize variables.
 */
static void ClearMeshIO(MeshIO *m) {
  m->io->pointlist                  = (REAL *)NULL;
  m->io->pointattributelist         = (REAL *)NULL;
  m->io->pointmarkerlist            = (int  *)NULL;
  m->io->edgelist                   = (int  *)NULL;
  m->io->edgemarkerlist             = (int  *)NULL;
  m->io->holelist                   = (REAL *)NULL;
  m->io->neighborlist               = (int  *)NULL;
  m->io->normlist                   = (REAL *)NULL;
  m->io->regionlist                 = (REAL *)NULL;
  m->io->segmentmarkerlist          = (int  *)NULL;
  m->io->segmentlist                = (int  *)NULL;
  m->io->trianglelist               = (int  *)NULL;
  m->io->triangleattributelist      = (REAL *)NULL;
  m->io->trianglearealist           = (REAL *)NULL;

  m->io->numberofpoints             = 0;
  m->io->numberofpointattributes    = 0;
  m->io->numberofcorners            = 0;
  m->io->numberofedges              = 0;
  m->io->numberofholes              = 0;
  m->io->numberofregions            = 0;
  m->io->numberofsegments           = 0;
  m->io->numberoftriangleattributes = 0;
  m->io->numberoftriangles          = 0;

  m->pointlistsize                  = 0; 
  m->pointattributelistsize         = 0;
  m->pointmarkerlistsize            = 0;
  m->edgelistsize                   = 0;
  m->edgemarkerlistsize             = 0;
  m->holelistsize                   = 0;
  m->neighborlistsize               = 0;
  m->normlistsize                   = 0;
  m->regionlistsize                 = 0;
  m->segmentmarkerlistsize          = 0;
  m->segmentlistsize                = 0;
  m->trianglelistsize               = 0;
  m->triangleattributelistsize      = 0;
  m->trianglearealistsize           = 0;
 }

/**
 * Create new instance.
 */
static PyObject* MeshIO_new(PyTypeObject* type, PyObject* args,
                            PyObject* kwds) {
  MeshIO* self;
  self = (MeshIO*)type->tp_alloc(type, 0);
  self->io = (struct triangulateio*) malloc(sizeof(struct triangulateio));
  ClearMeshIO(self);
  return (PyObject*)self;
}

/**
 * Initialize type.
 */
static int MeshIO_init(MeshIO *self, PyObject *args) {
  ClearMeshIO(self);
  return 0;
}

/**
 * Deallocate instance.
 */
static void MeshIO_dealloc(MeshIO* self) {
  free(self->io->pointlist);
  free(self->io->pointattributelist);
  free(self->io->pointmarkerlist);
  free(self->io->edgelist);
  free(self->io->edgemarkerlist);
  
  if(self->io->holelist != NULL)  {
    free(self->io->holelist);
    self->io->holelist = NULL;
  }
  
  free(self->io->neighborlist);
  free(self->io->normlist);
  
  if(self->io->regionlist!= NULL)  {
    free(self->io->regionlist);
    self->io->regionlist = NULL;
  }
  
  free(self->io->segmentmarkerlist);
  free(self->io->segmentlist);
  free(self->io->triangleattributelist);
  free(self->io->trianglearealist);
  free(self->io->trianglelist);
  free(self->io);
  self->ob_type->tp_free((PyObject*)self);
}


/*******************************************************************************
 * Functions for handling lists
 ******************************************************************************/
/**
 * Check if given PyList is ok.
 */
static int check_list(PyObject* list) {
  int ret = 0;
  /* Check if PyObject is NULL. */
  if (list == NULL) {
    PyErr_SetString(PyExc_TypeError, "error (1001): list pointer is NULL");
    ret = -1;
  }
  /* Check if it is indeed a list. */
  else if (!PyList_Check(list)) {
    PyErr_SetString(PyExc_TypeError, "error (1002): object is not a list");
    ret = -1;
  }
  return ret;
}

int set_double_list(PyObject* list, REAL** array, int* size,
                    int* numberofitems) {
  /* Check the list. */
  if (check_list(list) < 0) {
    return -1;
  }
  /* Create a new array, read data and check for consistency. */
  int n = PyList_Size(list);
  REAL* newarray = (REAL*)malloc(n*sizeof(REAL));
  int i;
  for (i = 0; i < n; i++) {
    PyObject* value = PyList_GetItem(list, i);
    if (PyFloat_Check(value)) {
      newarray[i] = PyFloat_AsDouble(value);
    } else {
      PyErr_SetString(PyExc_TypeError, "error (1011): object is not a float");
      free(newarray);
      return -1;
    }
  }
  /* Reading successful, assign values and return. */
  free(*array);
  *array = newarray;
  *size   = n;
  if (numberofitems != NULL) {
    *numberofitems  = 0;
  }
  return 0;
}

static PyObject* get_double_list(int size, REAL* array) {
  PyObject* list = PyList_New(size);
  int i;
  for (i=0; i<size; i++) {
    PyList_SetItem(list, i, Py_BuildValue("f", array[i]));
  }
  Py_INCREF(list);
  return list;
}

int set_int_list(PyObject* list, int** array, int* size,
                    int* numberofitems) {
  /* Check the list. */
  if (check_list(list) < 0) {
    return -1;
  }
  /* Create a new array, read data and check for consistency. */
  int n = PyList_Size(list);
  int* newarray = (int*)malloc(n*sizeof(int));
  int i;
  for (i = 0; i < n; i++) {
    PyObject* value = PyList_GetItem(list, i);
    if (PyInt_Check(value)) {
      newarray[i] = PyInt_AsLong(value);
    } else {
      PyErr_SetString(PyExc_TypeError, "error (1011): object is not an integer");
      free(newarray);
      return -1;
    }
  }
  /* Reading successful, assign values and return. */
  free(*array);
  *array = newarray;
  *size   = n;
  if (numberofitems != NULL) {
    *numberofitems  = 0;
  }
  return 0;
}

static PyObject* get_int_list(int size, int* array) {
  PyObject* list = PyList_New(size);
  int i;
  for (i=0; i<size; i++) {
    PyList_SetItem(list, i, Py_BuildValue("i", array[i]));
  }
  Py_INCREF(list);
  return list;
}

/*******************************************************************************
 * Points 
 ******************************************************************************/
static int MeshIO_setpoints(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->pointlist,
                               &self->pointlistsize,
                               &self->io->numberofpoints);
}
static PyObject* MeshIO_getpoints(MeshIO* self, void* closure) {
  return get_double_list(self->pointlistsize,
                         self->io->pointlist);
}
static PyObject* MeshIO_getnumberofpoints(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofpoints);
}


/*******************************************************************************
* Point attributes
 ******************************************************************************/
static int MeshIO_setpointattributes(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->pointattributelist,
                               &self->pointattributelistsize,
                               &self->io->numberofpointattributes);
}

static PyObject* MeshIO_getpointattributes(MeshIO* self, void* closure) {
  return get_double_list(self->pointattributelistsize,
                         self->io->pointattributelist);
}
static PyObject* MeshIO_getnumberofpointattributes(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofpointattributes);
}

/*******************************************************************************
 * Point markers 
 ******************************************************************************/
static int MeshIO_setpointmarkers(MeshIO* self, PyObject* list, void* closure) {
  return set_int_list( list, &self->io->pointmarkerlist,
                             &self->pointmarkerlistsize,
                             NULL);
}
static PyObject* MeshIO_getpointmarkers(MeshIO* self, void* closure) {
  return get_int_list(self->pointmarkerlistsize,
                      self->io->pointmarkerlist);
}

/*******************************************************************************
 * Triangles
 ******************************************************************************/
static int MeshIO_settriangles(MeshIO* self, PyObject* list, void* closure) {
  return set_int_list(list, &self->io->trianglelist,
                            &self->trianglelistsize,
                            &self->io->numberoftriangles);
}
static PyObject* MeshIO_gettriangles(MeshIO* self, void* closure) {
  return get_int_list(self->trianglelistsize,
                      self->io->trianglelist);
}
static PyObject* MeshIO_getnumberoftriangles(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberoftriangles);
}

/*******************************************************************************
 * Triangle attributes
 ******************************************************************************/
static int MeshIO_settriangleattributes(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->triangleattributelist,
                               &self->triangleattributelistsize,
                               &self->io->numberoftriangleattributes);
}
static PyObject* MeshIO_gettriangleattributes(MeshIO* self, void* closure) {
  return get_double_list(self->triangleattributelistsize,
                         self->io->triangleattributelist);
}
static PyObject* MeshIO_getnumberoftriangleattributes(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberoftriangleattributes);
}

/*******************************************************************************
 * Triangle areas
 ******************************************************************************/
static int MeshIO_settriangleareas(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->trianglearealist,
                               &self->trianglearealistsize,
                               NULL);
}
static PyObject* MeshIO_gettriangleareas(MeshIO* self, void* closure) {
  return get_double_list(self->trianglearealistsize,
                         self->io->trianglearealist);
}

/*******************************************************************************
 * Neighbors
 ******************************************************************************/
static PyObject* MeshIO_getneighbors(MeshIO* self, void* closure) {
  return get_int_list(self->neighborlistsize,
                      self->io->neighborlist);
}

/*******************************************************************************
 * Triangle corners 
 ******************************************************************************/
static PyObject* MeshIO_getnumberofcorners(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofcorners);
}

/*******************************************************************************
 * Segments
 ******************************************************************************/
static int MeshIO_setsegments(MeshIO* self, PyObject* list, void* closure) {
  return set_int_list(list, &self->io->segmentlist,
                            &self->segmentlistsize,
                            &self->io->numberofsegments);
}
static PyObject* MeshIO_getsegments(MeshIO* self, void* closure) {
  return get_int_list(self->segmentlistsize,
                      self->io->segmentlist);
}
static PyObject* MeshIO_getnumberofsegments(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofsegments);
}

/*******************************************************************************
 * Segment markers
 ******************************************************************************/
static int MeshIO_setsegmentmarkers(MeshIO* self, PyObject* list, void* closure) {
  return set_int_list(list, &self->io->segmentmarkerlist,
                            &self->segmentmarkerlistsize,
                            NULL);
}
static PyObject* MeshIO_getsegmentmarkers(MeshIO* self, void* closure) {
  return get_int_list(self->segmentmarkerlistsize,
                      self->io->segmentmarkerlist);
}

/*******************************************************************************
 * Holes
 ******************************************************************************/
static int MeshIO_setholes(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->holelist,
                               &self->holelistsize,
                               &self->io->numberofholes);
}
static PyObject* MeshIO_getholes(MeshIO* self, void* closure) {
  return get_double_list(self->holelistsize,
                         self->io->holelist);
}
static PyObject* MeshIO_getnumberofholes(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofholes);
}

/*******************************************************************************
 * Regions
 ******************************************************************************/
static int MeshIO_setregions(MeshIO* self, PyObject* list, void* closure) {
  return set_double_list(list, &self->io->regionlist,
                               &self->regionlistsize,
                               &self->io->numberofregions);
}

static PyObject* MeshIO_getregions(MeshIO* self, void* closure) {
  return get_double_list(self->regionlistsize,
                         self->io->regionlist);
}

static PyObject* MeshIO_getnumberofregions(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofregions);
}

/*******************************************************************************
 * Edges
 ******************************************************************************/
static PyObject* MeshIO_getedges(MeshIO* self, void* closure) {
  return get_int_list(self->edgelistsize,
                      self->io->edgelist);
}
static PyObject* MeshIO_getnumberofedges(MeshIO* self, void* closure) {
  return Py_BuildValue("i", self->io->numberofedges);
}

/*******************************************************************************
 * Edge markers
 ******************************************************************************/
static PyObject* MeshIO_getedgemarkers(MeshIO* self, void* closure) {
  return get_int_list(self->edgemarkerlistsize,
                      self->io->edgemarkerlist);
}

/*******************************************************************************
 * Norms
 ******************************************************************************/
static PyObject* MeshIO_getnorms(MeshIO* self, void* closure) {
  return get_double_list(self->normlistsize,
                      self->io->normlist);
}

/*******************************************************************************
 * Getters and setters
 ******************************************************************************/
static PyGetSetDef MeshIO_getseters[] = {
  {
    "points",
    (getter)MeshIO_getpoints,
    (setter)MeshIO_setpoints,
    "An array of point coordinates. The first point's x "
    "coordinate is at index [0] and its y coordinate at index [1], followed "
    "by the coordinates of the remaining points. Each point occupies two "
    "floats.",
    NULL
  },
  {
   "numberofpoints",
    (getter)MeshIO_getnumberofpoints,
    NULL,
    "Number of points",
    NULL
  },
  {
    "pointattributes",
    (getter)MeshIO_getpointattributes,
    (setter)MeshIO_setpointattributes,
    "An array of point attributes. Each point's "
    "attributes occupy `numberofpointattributes' floats.",
    NULL
  },
  {
    "numberofpointattributes",
    (getter)MeshIO_getnumberofpointattributes,
    NULL,
    "Number of point attributes.",
    NULL
  },
  {
    "pointmarkers",
    (getter)MeshIO_getpointmarkers,
    (setter)MeshIO_setpointmarkers,
    "An array of point markers; one int per point.",
    NULL
  },
  {
    "triangles",
    (getter)MeshIO_gettriangles,
    (setter)MeshIO_settriangles,
    "An array of triangle corners. The first triangle's "
    "first corner is at index [0], followed by its other two corners in "
    "counterclockwise order, followed by any other nodes if the triangle "
    "represents a nonlinear element. Each triangle occupies "
    "`numberofcorners' ints.",
    NULL
  },
  {
   "numberoftriangles",
    (getter)MeshIO_getnumberoftriangles,
    NULL,
    "Number of triangles",
    NULL
  },
  {
    "triangleattributes",
    (getter)MeshIO_gettriangleattributes,
    (setter)MeshIO_settriangleattributes,
    "An array of triangle attributes.  Each "
    "triangle's attributes occupy `numberoftriangleattributes' floats.",
    NULL
  },
  {
   "numberoftriangleattributes",
    (getter)MeshIO_getnumberoftriangleattributes,
    NULL,
    "Number of triangles attributes.",
    NULL
  },
  {
    "triangleareas",
    (getter)MeshIO_gettriangleareas,
    (setter)MeshIO_settriangleareas,
    "An array of triangle area constraints; one float per"
    "triangle. Input only.",
    NULL
  },
  {
    "neighbors",
    (getter)MeshIO_getneighbors,
    NULL,
    "An array of triangle neighbors; three ints per "
    "triangle. Output only.",
    NULL
  },
  {
   "numberofcorners",
    (getter)MeshIO_getnumberofcorners,
    NULL,
    "Number of triangle corners.",
    NULL
  },
  {
    "segments",
    (getter)MeshIO_getsegments,
    (setter)MeshIO_setsegments,
    "An array of segment endpoints. The first segment's "
    "endpoints are at indices [0] and [1], followed by the remaining "
    "segments. Two ints per segment.",
    NULL
  },
  {
   "numberofsegments",
    (getter)MeshIO_getnumberofsegments,
    NULL,
    "Number of segments",
    NULL
  },
  {
    "segmentmarkers",
    (getter)MeshIO_getsegmentmarkers,
    (setter)MeshIO_setsegmentmarkers,
    "An array of segment markers; one int per segment.",
    NULL
  },  
  {
    "holes",
    (getter)MeshIO_getholes,
    (setter)MeshIO_setholes,
    "An array of holes.  The first hole's x and y coordinates "
    "are at indices [0] and [1], followed by the remaining holes. Two "
    "REALs per hole.  Input only, although the pointer is copied to the "
    "output structure for your convenience.",
    NULL
  },
  {
   "numberofholes",
    (getter)MeshIO_getnumberofholes,
    NULL,
    "Number of holes",
    NULL
  },  
  {
    "regions",
    (getter)MeshIO_getregions,
    (setter)MeshIO_setregions,
    "An array of regional attributes and area constraints. "
    "The first constraint's x and y coordinates are at indices [0] and [1], "
    "followed by the regional attribute at index [2], followed by the "
    "maximum area at index [3], followed by the remaining area constraints. "
    "Four floats per area constraint. Note that each regional attribute is "
    "used only if you select the `A' switch, and each area constraint is "
    "used only if you select the `a' switch (with no number following), but "
    "omitting one of these switches does not change the memory layout. "
    "Input only, although the pointer is copied to the output structure for"
    "your convenience.",
    NULL
  },
  {
   "numberofregions",
    (getter)MeshIO_getnumberofregions,
    NULL,
    "Number of regions",
    NULL
  },
  {
    "edges",
    (getter)MeshIO_getedges,
    NULL,
    "An array of edge endpoints. The first edge's endpoints are "
    "at indices [0] and [1], followed by the remaining edges. Two ints per "
    "edge.  Output only.", 
     NULL
  },
  {
   "numberofedges",
    (getter)MeshIO_getnumberofedges,
    NULL,
    "Number of edges",
    NULL
  },
  {
    "edgemarkers",
    (getter)MeshIO_getedgemarkers,
    NULL,
    "An array of edge markers; one int per edge. Output"
    "only.",
    NULL
  },
  {
    "norms",
    (getter)MeshIO_getnorms,
    NULL,
    "An array of normal vectors, used for infinite rays in "
    "Voronoi diagrams. The first normal vector's x and y magnitudes are "
    "at indices [0] and [1], followed by the remaining vectors.  For each "
    "finite edge in a Voronoi diagram, the normal vector written is the "
    "zero vector. Two floats per edge. Output only.",
    NULL
  },
  {NULL, NULL, NULL, NULL}  /* Sentinel */
};

/**
 * Forward declaration. Needed in order PyObject_TypeCheck (see below) to work.
 */
static PyTypeObject trigen_MeshIOType;



static PyObject* MeshIO_mesh(PyObject *self, PyObject *args) {
  char *switches;
  PyObject* min = NULL;
  PyObject* mout = NULL;
  if (!PyArg_ParseTuple(args, "sOO", &switches, &min, &mout)) {
    PyErr_SetString(PyExc_TypeError, "Expected (string, MeshIO)");
    return NULL;
  }
  if (!PyObject_TypeCheck(min, &trigen_MeshIOType)) {
    PyErr_SetString(PyExc_TypeError, "Argument 2 is not a Meshio");
    return NULL;
  }

  MeshIO* in = (MeshIO*)min;

  /* number of points */
  if (in->pointlistsize%2 !=0) {
    PyErr_SetString(PyExc_TypeError, "error (1031): point error.");
    return NULL;
  } else {
    in->io->numberofpoints = in->pointlistsize/2;
  }

  /* number of point attributes */
  if (in->pointattributelistsize%(in->pointlistsize/2) !=0) {
    PyErr_SetString(PyExc_TypeError, "error (1032): point attribute error.");
    return NULL;
  } else {
    in->io->numberofpointattibutes = 
      in->pointattributelistsize/(in->pointlistsize/2);
  }

  /* number of corners */
  /* TODO:not implemented yet*/

  /* number of edges */
  /* TODO:not implemented yet*/

  /* number of triangles */
  /* TODO:not implemented yet*/
  
  /* number of triangle attributes */
  /* TODO:not implemented yet*/
  
  /* number of segments*/
  if (in->segmentlistsize%2 !=0) {
    PyErr_SetString(PyExc_TypeError, "error (1037): segment error.");
    return NULL;
  } else {
    in->io->numberofsegments = in->segmentlistsize/2;
  }

  /* number of holes */
  if (in->holelistsize%2 !=0) {
    PyErr_SetString(PyExc_TypeError, "error (1038): hole error.");
    return NULL;
  } else {
    in->io->numberofholes = in->holelistsize/2;
  }

  /* number of regions */
  if (in->regionlistsize%4 !=0) {
    PyErr_SetString(PyExc_TypeError, "error (1039): region error.");
    return NULL;
  } else {
    in->io->numberofregions = in->regionlistsize/4;
  }
  
  MeshIO* out = (MeshIO*)mout;
 
  triangulate(switches, in->io, out->io, (struct triangulateio*)NULL);

  
  out->pointlistsize   = 2*(out->io->numberofpoints);
  out->segmentlistsize = 2*(out->io->numberofsegments);

  out->holelistsize    = 2*(out->io->numberofholes);
  out->io->holelist    = (REAL*) malloc(out->holelistsize * sizeof(REAL));

  int i;
  for (i = 0; i < out->holelistsize; i++) {
    out->io->holelist[i] = in->io->holelist[i];
  }

  out->regionlistsize  = 4*(out->io->numberofregions);
  out->io->regionlist  = (REAL*) malloc(out->regionlistsize * sizeof(REAL));
  for (i = 0; i < out->regionlistsize; i++) {
    out->io->regionlist[i] = in->io->regionlist[i];
  }
  out->triangleattributelistsize = (out->io->numberoftriangles)*(out->io->numberoftriangleattributes);
  return Py_None;
}  







static PyObject* MeshIO_vtk(PyObject *self, PyObject *args) {
  char *filename;
  PyObject* mio = NULL;
  if (!PyArg_ParseTuple(args, "sO", &filename, &mio)) {
    PyErr_SetString(PyExc_TypeError, "Expected (string, MeshIO)");
    return NULL;
  }
  
  MeshIO* m = (MeshIO*)mio;
  FILE* fp = fopen(filename, "w");
  fprintf(fp, "# vtk DataFile Version 3.0\n");
  fprintf(fp, "nemesis vtk output\n");
  fprintf(fp, "ASCII\n");
  fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
  fprintf(fp, "\n");
  fprintf(fp, "POINTS %d float\n", m->io->numberofpoints);
  int i;
  for (i=0; i<m->io->numberofpoints; i++) {
    double x = m->io->pointlist[2*i];
    double y = m->io->pointlist[2*i+1];
    fprintf(fp, "%.4f %.4f %.4f\n", x, y, 0.);
  }
  fprintf(fp, "\n");
  fprintf(fp, "CELLS %d %d\n", m->io->numberoftriangles, 4*(m->io->numberoftriangles));
  for (i=0; i<m->io->numberoftriangles; i++) {
    int n1 = m->io->trianglelist[3*i]  -1;
    int n2 = m->io->trianglelist[3*i+1]-1;
    int n3 = m->io->trianglelist[3*i+2]-1;
    fprintf(fp, "3 %d %d %d\n", n1, n2, n3);
  }
  fprintf(fp, "\n");
  fprintf(fp, "CELL_TYPES %d\n", m->io->numberoftriangles);
  for (i=0; i<m->io->numberoftriangles; i++) {
    fprintf(fp, "5\n");
  }

  fclose(fp);
  return Py_None;
}

static PyMethodDef trigen_methods[] = {
  {"mesh", (PyCFunction)MeshIO_mesh, METH_VARARGS, "Mesh."},
  {"vtk",  (PyCFunction)MeshIO_vtk,  METH_VARARGS, "Export to vtk output format."},
  {NULL, NULL, 0, NULL}  /* Sentinel */
};

static PyTypeObject trigen_MeshIOType = {
PyObject_HEAD_INIT(NULL)
  0,                          /* ob_size            */
  "trigen.MeshIO",            /* tp_name            */
  sizeof(MeshIO),             /* tp_basicsize       */
  0,                          /* tp_itemsize        */
  (destructor)MeshIO_dealloc, /* tp_dealloc         */
  0,                          /* tp_print           */
  0,                          /* tp_getattr         */
  0,                          /* tp_setattr         */
  0,                          /* tp_compare         */
  0,                          /* tp_repr            */
  0,                          /* tp_as_number       */
  0,                          /* tp_as_sequence     */
  0,                          /* tp_as_mapping      */
  0,                          /* tp_hash            */
  0,                          /* tp_call            */
  0,                          /* tp_str             */
  0,                          /* tp_getattro        */
  0,                          /* tp_setattro        */
  0,                          /* tp_as_buffer       */
  Py_TPFLAGS_DEFAULT,         /* tp_flags           */
  "MeshIO object",            /* tp_doc             */
  0,                          /* tp_traverse        */
  0,                          /* tp_clear           */
  0,                          /* tp_richcompare     */
  0,                          /* tp_weaklistoffset  */
  0,                          /* tp_iter            */
  0,                          /* tp_iternext        */
  trigen_methods,             /* tp_methods         */
  0,                          /* tp_members         */
  MeshIO_getseters,           /* tp_getset          */
  0,                          /* tp_base            */
  0,                          /* tp_dict            */
  0,                          /* tp_descr_get       */
  0,                          /* tp_descr_set       */
  0,                          /* tp_dictoffset      */
  (initproc)MeshIO_init,      /* tp_init            */
  0,                          /* tp_alloc           */
  MeshIO_new,                 /* tp_new             */
};

/* declarations for DLL import/export */
#ifndef PyMODINIT_FUNC
#define PyMODINIT_FUNC void
#endif

PyMODINIT_FUNC inittrigen(void) {
  PyObject* m;
  if (PyType_Ready(&trigen_MeshIOType) < 0) {
    return;
  }
  m = Py_InitModule3("trigen", trigen_methods,
                     "Wraper module for the triangle library.");
  Py_INCREF(&trigen_MeshIOType);
  PyModule_AddObject(m, "MeshIO", (PyObject*)&trigen_MeshIOType);
}

